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19.  Abstract 

The  research  supported  by  this  grant  was  essentially  concerned  with  the  numerical  solution  of 
the  Zakai  equation  from  nonlinear  filtering.  The  Zakai  equation  is  a  linear  stochastic  equation  of 
parabolic  type  and  its  solution  provides  a  nonnormalized  probability  density  from  which  we  can 
reconstruct  the  solution  of  a  dynamical  system  governed  by  stochastic  ordinary  differential  equations 
when  the  observation  is  itself  noisy.  The  Zakai  equation  presents  some  unique  features  in  the  sense 
that 

(i)  It  has  practical  applications; 

(ii)  It  is  a  highly  advective,  advention-diffusion  linear  parabolic  equation; 

(iii)  The  space  dimension  is  very  large  (of  the  order  of  ten  and  more  for  practical  applications); 

(iv)  It  is  a  stochastic  equation. 

Each  of  the  difficulties  listed  in  (ii),  (iii),  (iv)  leads  to  very  challenging  problems  when  it  goes  to  the 
numerical  solution,  but  their  combination  provides  a  formidable  numerical  problem  which  from  our 
point  of  view  requires  computer  power  beyond  the  possibility  of  the  existing  hardwares  and  softwares. 

In  view  of  developing  know  now  we  have  been  investigating  systematicc  methods  for  solving  the  Zakai 
equation,  which  have  proved  reliable  for  solving  problems  with  space  dimension  six  at  most  (this 
number  will  increase  with  the  progress  of  computers).  The  methodology  that  we  have  been  studying 
relies  on  operator  splitting  methods  in  which  one  decouples,  via  time  discretization,  the  stochastic  part 
and  the  deterministic  part  of  the  Zakai  equation.  Using  this  approach  we  obtain  subproblems  which 
have  either  closed  form  solution  (the  stochastic  part)  or  can  be  solved  (the  deterministic  part;  modulo 
the  difficulty  of  the  dimension)  by  efficient  methods  for  traditional  elliptic  and  parabolic  problems. 
The  same  methodology  can  be  extended  to  the  parabolic  inequalities  of  obstacle  type  originating  from 
impulse  control. 

Concerning  the  space  approximation  itself  it  can  be  achieved  by  finite  difference  or  finite  element 
methods,  with  a  specific  treatment  of  the  first  order  terms  for  which  the  first  order  upwinding  methods 
i  investigated  by  some  authors  are  by  far  too  dissipative.  We  have  been  therefore  investigating  second 

i 

order  upwinding  methods  which  are  by  far  much  less  dissipative  (in  the  detailed  final  report  we  also 
comment  on  third  order  upwinding  methods). 
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Title:  Numerical  Methods  for  Parabolic  Equations  ic  Inequalities  in  Very  High  Dimensions. 


1.  Introduction.  Motivation 

)  The  main  goal  of  this  projects  was  to  investigate  the  numerical  solution  of  parabolic  equations 
in  very  high  dimension ,  the  principal  motivation  being  the  solution  of  the  so-called  Zakai  equation  from 
nonlinear  filtering.  As  it  will  be  seen  in  Section  2,  below,  this  equation  is  of  the  adveciion- diffusion 
type  and  therefore  is  closely  related  to  other  equations  of  this  type  such  as  the  Navier-Stokes  equations 
from  Fluid  Dynamics;  since  these  equations  are  also  of  interest  for  the  Department  of  Defense  some 
aspects  of  their  numerical  solution,  byproducts  of  the  present  investigation  will  be  considered  in  the 
present  document. 

The  content  of  this  report  is  the  following: 

In  Section  2,  we  describe  the  Zakai  equation  and  its  origin.  In  Section  3,  we  consider  the  solution  of 
time  dependent  problems  by  operator  splitting  methods,  and  show  in  Sections  4,  5,  6  how  these 
methods  apply  to  the  Zakai  equation,  to  parabolic  inequalities  of  obstacle  type  and  to  Navier-Stokes 
equations.  N 


‘In  Section  7  we  comment  on  the  solution  of  advection-diffusion  problems  by  upwinding,  particle  and 
characteristics  methods.  In  Section  8  we  go  back  to  the  Zakai  equation  and  its  stochastic  dynamical 
system  origin  and  comment  about  the  feasibility  of  the  Zakai  equation  approach  for  the  analysis  of 
such  systems.  / 


I  v K  )  C~ 


2.  The  Zakai  equation. 

We  consider  a  signal  process  given  by 
(2.1)  dx,  =  b(t,  x,)dt  -(-  cr(t,  x<)dw,. 
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Here  t6[0,T],  xERd  and  x0  has  a  given  probability  density  p0(x).  The  vector  w,  =  |w  1 . wf} is  a 

d-dimensional  Brownian  motion  defined  on  a  probability  space  (Q,  F,  p*).  Also,  b(t,  x)  = 

{bj(t,  x)j  and  <r(t,  x)=  (<r,-y(t,x)),  i J  =  1,  d,  have  components  defined  on  [0,T]xRd  which  are 

bounded  and  smooth.  We  suppose  now  that  the  observation  is  given  by 

(2.2)  dy,=  h(t,  x()  +  g(t)dwt  +  g(t)dwt. 

Here,  y€Rm,  yo  =  0  and  wt  =  jwj  j  is  a  m  dimensional  Brownian  motion  defined  on  (Q,  F,  p*)  and 

is  independent  of  w.  The  function  h(t,  x)  =  jhy(t,  x)|^  has  essentially  the  same  smoothness 
and  boundedness  properties  than  b  and  cr.  We  shall  suppose  that  matrices  g(t)  =  (g.  (t)), 
i=l,  ...  m,  j=l,  ...  d,  and  g(t)  =  (gt,(t)),  k,  1=1,  ...m  are  continuous  on  [0,T].  We  suppose  that 

(2-3)  g(t)g(t)T+g(t)g(t)T=  I. 

Condition  (2.3)  is  a  normalization  hypothesis  which  can  always  be  satisfied  by  re-scaling  the 
observation  process. 


Denoting  by  g y  the  j-th  row  of  g  we  define  the  vector  c;-  by 
C>(t)  =  <r(t)  gT(t), 

and  then  h-  by 

h;(t)  =  hy(t)  -  div  Cy(t),  j=l,  ...  d. 


Define  now  an  operator  Hy,  mapping  H^R®)  into  L2(Rd),  by 


Hy(t)  u(t)  =  hy(t)  u(t)  -  (Cy(t),  Vu(t)), 


for  u  G  H^R"). 


The  elliptic  operator  L*  is  the  formal  adjoint  of  the  generator  of  (2.1),  i.e. 


here 


-  {.,}  = 


T 

aa  . 
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Define  Y,  =  cr|y,:  0  <  s  <  t|  for  the  cr-field  generated  by  the  observations  up  to  time  t. 
There  exists  an  unnormalized  density  q,  of  xt  conditioned  on  Y(;  q,  satisfies 

(2.4)  qQ  =  pQ 

and  the  following  stochastic  partial  differential  equation 

t  t 

(2.5)  q,(x)  =  Pq(x)  +  [  L*  q,(x)  ds  +  £  [  Hj(s)  Q»(x)  d4- 

0  j=1  0 

for  t  £  [0,T].  Equation  (2.5)  is  precisely  the  Zakai  equation  (see,  e.g.,  [5],  [10]  for  more  details). 


3.  Operator  Splitting  Methods  for  Time  Dependent  Problems. 

Let’s  consider  a  (deterministic)  dynamical  system  modelled  by 

(3.1)  ^  +  A(u)  =  0,  for  t>0, 

(3.2)  u(0)  =u0. 

In  (3.1),  (3.2),  A  is  a(possibl>  nonlinear)  operator  from  a  Hilbert  space  H  into  itself,  and  u0eH.  We 
suppose  that  A  has  a  nontrivial  decomposition  of  the  following  type 

(3.3)  A  =  Aa  +  A2; 

by  nontrivial  we  mean  that  Aa,  A2  are  individually  simpler  than  A  (indeed  Aa  and/or  A2  can  be 
multivalued  operators). 

There  exist  many  schemes  taking  advantage  of  decomposition  (3.3)  (see,  e.g.,  [1],  [2],  [3]  for  such 
schemes).  We  shall  present  here  some  of  them: 
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3.1  Fractional  Step  Schemes 


1st.  Scheme: 


(3.4) 


u  =  u0; 


then  for  n>0,  being  known  we  compute  u"+1^2  and  then  un+I  as  follows 


(3.5) 


H"Tt~  +  A1(u,l+1/2)  =  0, 


„n  +  l  n+ 1/2 

(3.6)  U - =jJ! -  +  A2(un)  =  0. 


This  scheme  is  at  most  first  order  accurate  in  At.  To  obtain  a  better  accuracy,  we  can  “symmetrize” 
scheme  (3.4)  -  (3.6)  by  using  the  following  scheme 

2nd  Scheme: 


(3'7)  +  Al  (un+1/3)  =  °’ 


(3.8) 


n+2/3  n+l/3 

U  —  U 


+  A2  (un+2/3)  =  0, 


At 


,,n  +  l  n  +  2/3 

(3-9)  ~ — Xl72 -  +  Ai  <u  >  =  °- 


It  follows  from,  e.g.,  [3]  (see  also  [4]),  that  in  many  situations,  the  above  scheme  is  second  order 
accurate  with  respect  to  At. 


3.2.  Alternating  Direction  Schemes. 

The  above  schemes  are  not  too  well  suited  to  capture  steady  state  solutions  of  (3.1),  i.e., 
solutions  of 
(3.10)  A(u)  =  0 
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(if  such  solutions  exist).  This  drawback  is  corrected  if  one  uses  the  schemes  below: 


Peaceman-Rachford  Scheme: 

(3.10)  u°  =  u0; 

then  for  n  >  0, 

(3-H)  "n+At/f-n  +  Al  (U"+1/2)  +  A^u")  =  °> 

(3-12)  +  Al  (U”+1/2)  +  A2  (U"+1)  =  °‘ 

Douglas-Rachford  Scheme: 

(3.13)  u°  =  u0; 
then  for  n  >  0, 

(3.14)  y--2f  +  Ai  (U"41/2)  +  A2(u")  =  °- 

(3.15)  u"^--un  +  A,  (u"+1/2)  +  A2  (u°+1)  =  0. 

fl-Scheme:  With  9  €  (0,,5),  this  scheme  is  defined  by  (3.13)  and 
(•U6)  u-By-t-u"  4-  A,  (u"+(?)  +  A2(un)  =  0, 

,,n  +  l— 0  n  +  6  n  .a 

(3‘17)  (l~29j  +  Al  (U  }  +  A*<u  >  =  °’ 

(3.14)  +  Al  („«+»)  +  A2(un+1-0)  =  0. 
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All  these  schemes  are  first  order  accurate  with  respect  to  At  in  general,  the  last  one  eing  (for  0  well 
chosen)  better  suited  than  the  first  two  to  capture  steady  states. 


4.  Application  to  the  solution  of  the  Zalcai  equation. 

Using  the  notation  of  Section  2,  we  are  considering  the  solution  of 

(4.1)  dy  4-  A*(t)y  dt  =  B(t)  y  •  dw. 

The  idea  here,  is  to  consider 

(4.2)  A*(t)y  dt  -  B(t)y  •  dw 

as  the  sum  of  two  operators  and  using  the  approach  associated  to  scheme  (3.4)  -  (3.6),  consider  a 
sequence  of  problems  of  the  form 

(4.3)  dy>  +  A(t )ip  dt  =  0, 

(4.4)  d\k  =  B(t)^  •  dw, 

which  are  considerably  simpler  than  (4.1);  indeed,  the  ip  equation  is  deterministic  and  the  9  equation 
has  a  closed  form  solution. 

The  algorithm  investigated  has  been  analyzed  in  [5]  and  generalized  to  more  complicated  stochastic 
partial  differential  equations  in  [6];  it  is  described  as  follows: 

Consider  first  (4.1)  on  the  time  interval  [0,T]  and  define  At=  k  -  T/(N+1).  We  shall  define  two 
processes  yu,  y2i  depending  on  k.  We  split  therefore  [0,  T]  in  [0,  k[,  [k,  2k[,  ...  [Nk,  (N+i)k[. 
Consider  now  [rk,  (r+l)k[  with  r=0,  1,  ...,  N;  then  yu,  y2k  are  defined  on  this  interval  by  the  relations 

(4.5)  dylt  +  A*(t)  yljfcdt  =  0, 
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(4.6) 


dy2i  =  B  (i)  yn  •  dw, 


(4.7)  yu(rk)  =  y'; 

(4.8)  y2t(rk)  =  y/+1/2 

the  sequences  yrk ,  y  k+l/'2  are  defined  as  follows 

(4.9)  y  k+1'2  =  yu  ((r+l)  k-o), 

(4.10)  yir+1  =  y2fc  ((r+1)  k-0). 

Relations  (4.5)  -  (4.10)  define  completely  ylk,  y2k  on  [0,  T[.  The  processes  ylk,  y2k  are  right 
continuous  and  their  discontinuity  points  are  k,  ...  Nk  (on  [0,  T[);  we  observe  that  (4.5)  is 
deterministic. 


The  convergence  of  yu,  y2i  to  the  solution  y  of  (4.1)  is  proved  in  [5],  justifying  therefore  the 
splitting  approach. 

In  practice,  problem  (4.5)  is  an  advection-diffusion  problem  whose  solution  will  be  discussed  in 
Section  7. 


5.  Application  to  the  solution  of  obstacle  type  parabolic  inequalities. 

It  follows  from  [7],  that  applications  in  impulse  control,  lead  to  the  solution  of  parabolic 
inequalities  of  the  following  type 

(5.1)  +  Au  -  f)  (tf  -  u)  =  0,  a.  e.  on  nx(0,T), 

(5.2)  +  Au  -  f  <  0,  a.  e., 

(5.3)  u  —  <  0,  a.e., 
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(5.4)  u  =  u0  at  t  =  0; 

in  (5.1)  (5.4),  A  is  an  advection-diffusion  elliptic  operator,  and  defines  the  obstacle.  From  a 

numerical  point  of  view,  it  may  be  more  convenient  to  formulate  (5.1)  -  (5.4)  as  a  parabolic  variational 
inequality  such  as 

(5.5)  [  ^4-  +  Au  —  fj  (v  —  u)  dx  >0,  VveK(t), 

a 

(o.G)  u(t)  (E  K(t),  Vt  >  0, 

(5.7)  u(0)  =  u0  (e  K  (0)), 
with 

(5.8)  K(t)  =  {v|v€V,  v(x,t)  <  *  (x,t)  a.t.  on  nx(0,T)}; 
in  (5.8),  V  is  a  Hilbert  space  of  Sobolev  type. 


In  (5.5)  -  (5.8)  we  have  implicitelv  assumed  that  'P,  and  therefore  the  (convex)  set  K,  vary 
continuously  with  t.  Concerning  the  numerical  solution  of  (5.5)  -  (5.8)  by  operator  splitting  methods, 
let’s  introduce 

(5.9)  C(t)  =  |v|v£L2(fi),  v(x)  <  9/(x,t)  a.e.  on  fij, 
and  the  indicator  functional  Ic  of  C(t),  defined  by 

(5.10)  Ic(v,  t)  =  0  i/v  6  C(t),  Ic(v,t)  =  +oo  i/v  g  C(t). 

Problem  (5.5)  -  (5.8)  can  then  be  formulated  as  a  nonlinear  parabolic  uequaiionr'  (it  is  in  fact  a 
multivalued  equation)  by 

(5.11)  0  €  |^  +  Au  +  0Ic(u)  -  f, 

(5.12)  u(0)  =  u0; 
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in  (5.11),  0Ic(u)  is  the  subdifferential  of  the  convex  functional  Ic  at  u;  it  is  defined,  Vv£L2(Q),  by 

(5.13)  wedlc(v,t)  O  Ic(z,t)  —  Ic(v,t)  >  | (z-v)  w  dx,  Vz£L2(Q). 

q 

Problem  (5.11),  (5.12)  can  be  approximated  by  the  following  nonlinear  parabolic  equation  (with  e>0): 

(5.14)  +  Au£  +  l(u£-*)+  =  f, 

(5.15)  u£(0)  =  u0, 

which  corresponds  to  a  penalty  approximation  of  condition  (5.6).  Operator 

(5.16)  v  — ►  (v  -  \&)+ 

is  (for  t  given)  a  (“nice”)  monotone  operator  from  L2(Q)  ->  L2(Q). 


Both  problems  (5.11),  (5.12)  and  (5.14),  (5.15)  can  be  solved  by  operator  splitting  methods.  For 
simplicity,  we  shall  just  describe  the  application  of  scheme  (3.4)-(3.6),  but  the  other  schemes  described 
in  Section  3  also  apply.  We  obtain  then  for  problem  (5.11),  (5.12) 

(5.17)  u°  =  u0; 


then  for  n>0,  we  compute  u"  +  1/2  and  un+1  via 


n  +  1/2 

(5.18)  ^  A~  -1-  +  die  (un+1/2,  (n+l)At)  =  0, 


(5.19) 


At 

un+1-  u"*1/2 


n  + 1  rn  + 1 


At 


+  Au  =  f 


Problem  (5.19)  is  an  advection-diffusion  elliptic  problem  whose  solution  is  discussed  in  Section  7;  on 
the  other  hand  problem  (5.18)  has  a  closed  form  solution  which  is  given  by 

(5.20)  u"+1/2(x)  =  min  (u"(x),  ^(x,  (n+l)At)). 
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A  similar  conclusion  would  hold  for  problem  (5.14),  (5.15).  In  fact  from  the  simplicity  of  scheme 
(5.17),  (5.19)  it  seems  that  penalty  is  useless  in  this  context;  we  have  mentioned  it  however  since  it 
may  be  useful  in  some  other  situations.  We  can  switch  of  course  the  splitting  order  and  solve  first  the 
advection  diffusion  part,  and  then  project  as  in  (5.20). 


6.  Application  to  the  Navier-Stokes  equations. 

Another  advection-diffusion  problem  of  interest  is  associated  to  the  Navier- Strokes  equations 
for  viscous  fluids.  Concentrating  on  the  incompressible  case  we  obtain 

(6.1)  -  i/Au  +  (u  •  V)  u  +  Vp  =  f  in  Cl. 

(6.2)  V-u  =  0  in  Cl  ( incompressibility  condition), 

(6.3)  u  (x,  0)  =  u0  (x)  in  Cl, 

with  appropriate  boundary  conditions;  for  simplicity  we  shall  assume  that 

(6.4)  u  =  g  on  r(=3fl). 

In  (6.1)  -  (6.4),  u  =  {u,.},.^  is  the  velocity,  p  is  the  pressure,  Cl  is  the  flow  region,  i/(>0)  is  a  viscosity 
coefficient  and  f  a  density  of  external  forces  and 

(6.5)  (v-V)w  =  {£  v,  . 

Following  the  approach  in  [8]  (see  also  [9])  we  shall  discretize  (6.1)  -  (6.4)  by  the  following  8-scheme, 
which  decouples  nonlinearity  and  incompressibility : 

(6.6)  u°  =  u0, 
and  then  for  n>0, 

(6.7)  i  -  auAnn+e  +  Vpn+<?  =  f"+(?  +  Aun  -  (un-V)un  in  Cl, 
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(6.7) 2  V-un+^  =  0  in  ft, 

(6.7) 3  un+*  =  gn+9  on  r, 


(6.8),  -!!  itth—  -/?^Aun+1-^+(un+1-0-V)  un+1-^=r+1-^  +  al/Aun+0-Vpn+(?  in  ft, 


'1  (l-20)At 


(6.8)2 


n  +  1-0  n  +  l-$  r 

u  =  g  on  r, 


n+1  ,  n+1-6 

(6.9),  U - — 1 - 


0At 


-ai/Aun+1+Vpn+1=  r+14-/3i/Aun+1~0-(un+1~6'-V)un+1_f?  in  ft, 


(6.9)9  V.un+1  =  0  in  ft, 


(6.9)3  un+1  =  gn+1  on  T. 

A  good  choice  for  6  is  6  =  1-1/42  and  for  a  and  /? 

0  =  jL  a  =  L2£ 
a  1-0'  P  1-0 - 

Numerical  results  obtained  using  the  above  scheme,  combined  to  finite  element  methods  are  given  in 
[8];  the  above  methods  has  been  extended  to  two-phase  flow  problems  in  [9]. 


7.  Numerical  Solution  of  Advection  -  Diffusion  Problems 

7.1.  Synopsis.  Motivation 

We  follow  here  the  presentation  in  [8],  completed  by  some  recent  findings  concerning 
upwinding  methods  of  order  three.  In  order  to  concentrate  on  Zakai  equation  related  problems  we  shall 
illustrate  the  subsequent  developments  by  considering  an  advection-diffusion  problem  associated  to  the 
following  stochastic  ordinary  differential  equation 
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(7.1) 


dx  =  V(x)  dt  +  dw, 


where  x  is  an  N-dimensional  vector  and  w  a  noise.  Assuming  convenient  hypotheses  on  the  noise  w, 
the  probability  of  finding  at  time  t  the  state  vector  x,  in  neighborhood  of  of  measure  d£  = 

d^...d^N,  is  p(£,t)  d£  where 


(7.2) 


P(^t) 


q(£,t) 

q(rj,t)dr? 


where  the  unnormilizaed  probability  density  q  satisfies  a  parabolic  equation  (see  Section  2  for  details). 
In  the  particular  case  where  V  is  divergence  free  and  for  simple  noise  models,  this  parabolic  equation 
reduces  to 


(7.3)  ^9  —  e  V2q  +  V  •  Vp  =  0. 

An  interesting  and  difficult  case  is  the  one  where  the  level  of  noise  is  “weak”,  implying  that  €  is 
“small”.  We  have  then  an  advection  dominated  adveciion- diffusion  equation. 


Using  the  splitting  methods  described  in  Section  4,  the  solution  of  such  equations  plays  a  fundamental 
role  in  the  implementation  of  some  solution  methods  for  the  Zakai  equation. 


In  the  sequel,  we  consider  the  following  initial  boundary  value  problem. 

(7.4) j  _  ev2p  +  V  •  Vp  =  f  in  f2x(0,T), 

(7.4) 2  p  =  g  on  dtt  x  (0,T), 

(7.4) 3  p(x,0)  =  p0  (x)  in  fi, 

with  fl  C  RN. 

For  solving  such  problems,  particularly  the  ones  associated  to  filtering,  one  has  to  face  two  outstanding 
difficulties,  namely 
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(i)  When  £  is  small,  the  problem  is  advection- dominated, 

(ii)  For  practical  problems,  we  have  N>3  (in  fact  we  may  have  N~10  and  more). 

The  methods  described  here  have  been  applied  in  [8]  to  N=2  to  6. 

7.2  Solution  of  problem  (7.4)  by  finite  differences  and  upwinding  methods. 

The  methods  we  need  have  to  be  accurate  and  robust-,  in  our  opinion,  we  have  to  avoid  those 
methods  requiring  parameter  tuning. 


For  simplicity,  we  consider  the  case  where  fl  =  (0,1)^,  with  N=2;  the  extension  to  N>2  is 
straight  forward.  With  I  a  positive  integer ,  we  define  h  by  h=l/I+l  and  consider  over  Q=Q UdQ  the 
mesh  points 

(7.5)  Mu  =  {ih,  jh};  0  <  ij<I+l. 

At  the  points  M,-;-  interior  to  fi  (i.e.,  1  <i J <I)  we  approximate  (7.4)  by  the  following  finite  difference 
scheme  (with  V={V1,V2}): 


p"*1  -  p?i  p/+ii  +  p?-+u  +  p,"++i+  p,"-\  -  ^.r1 
At  h2 


(7-6)!  { 


,  pn+1-p.n1+1 

+  V|  (M,y)  -ii - £=U- 


Pn  +  1  -  p."-*-,1 

+  V+  (MtJ) 


pn+\  -  pn+1 

V"  (M0)  -  ^  P,-l  - 


n  +  l  n+1 

(M,7)  -  -**1  h  P<i 


=  (n+l)At), 


with  ,  in  (7.6),  At(>0)  a  time  discretization  step,  p ~  p(M,;-,  nAt),  at  =  max  (0,a),  a+  =  max 
(0,-a),  Va  €  R,  and 


(7.6)2  p^,+  1  =  g(M*„  (n+1)  At)  if  M*,  €  dtt, 
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pkl  =  Po^M*')' 


(7.G)3 

Scheme  (7.6)  is  of  the  backward  Euler  type  for  the  time  discretization  and  of  the  first  order  upwinded 
type  for  the  space  discretization.  Probabilists  favor  the  above  finite  difference  scheme  because  it 
satisfies  a  discrete  maximum  principle  and  therefore  possesses  a  probabilistic  interpretation. 
Unfortunately  the  above  scheme  is  only  first  order  accurate ,  very  dissipative  and  not  well-suited  for 
those  situations  where  t  is  small  and  V  has  a  fast  variation  over  fi. 

An  interesting  alternative  to  (7.6)  is  obtained  through  a  space/time  discretization  which  is 
second  order  accurate  and  also  of  the  upwind  type  (however,  it  does  not  satisfy  the  discrete  maximum 
principle).  Such  a  scheme  is  obtained  as  follows: 

(7.7)i  pkl  =  pkl  15  °bta*ned  ( f°T  example)  via  (7.6); 

then  for  n>l  and  2<ij<I-l,  discretize  (7.4) ^  by 


(7.7)2 


If  M(J-  E  ft  with  either  M^.^j  or  belonging  to  T,  it  is  possible  that  or  .  are 

outside  fl;  in  such  a  case  we  can  use  to  discretize  V  -Vp  at  M,y,  a  first  order  scheme  like  in  (7.6) j,  or 
alternatively  a  centered  second  order  approximation  like 

(7.8)  (V- Vp)  (M{j)  ~  V1(M,J)  --j2HP--J'  +  V2(M0.)  -  P-;+12~  Pii~l 


-  n"+1  -?n"  4-1  nn_1 
2  P  ij  Zp.j+2PU' 


1  p.  *1  +  1  .  n  +  1  ,  n  +  1  4  —  Ti  +  l 

P.'  +  lj  +  P.-lj  +  P.;  +  l  +  P.j-1  ~  4p,-; 


3  0"+!  _  On  n+1  4.  !nn  +  1  3  nn+1—  9nn+1  4-lnn+1 


3^0  +  1  o_.  n  +  1  ,  1  n  +  1  3^  n  +  1 

+  V+  <M„)  ■  IJ;‘  2Pi'-;  +  V-  <M0) 


9d’1  +  1  4.  iD"  +  l 
*tj+l  ^  2P'j  +  2 

h 


=  (n-f-l)At). 
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The  boundary  conditions  are  treated  as  in  (7.6)2-  The  fact  that  the  problems  under 
consideration  may  have  a  fast  dynamics,  requires  the  use  of  small  h  and  At;  indeed,  we  can  take 
advantage  of  the  fact  that  At  is  small,  to  solve  the  above  discrete  elliptic  problems  by  successive  over¬ 
relaxation,  since  the  method  has  good  vectorization  and  parallelization  properties  (in  practice  few 
iterations  will  insure  convergence  at  each  time  step). 

Numerical  experiments  (see,  e.g.,  [8])  show  the  superiority  of  the  second  order  upwinding 
method,  over  the  first  order  method  (it  is  more  accurate,  less  dissipative  and  almost  as  easy  to 
implement). 

7.3.  A  third  order  in  space  upwinded  scheme. 

Influenced  by  third  order  in  space,  upwinded  finite  difference  schemes,  recently  developed  in 
Japan,  for  Computational  Fluid  Dynamics  purposes,  we  shall  describe  here  a  variant  of  scheme  (7.7) 
which  uses  a  third  order  accurate  space  discretization  of  the  advection  and  which  is  second  order 
accurate  in  time.  The  basic  principle  of  this  scheme  is  to  combine  two  second  order  accurate  schemes, 
one  centered  (in  fact,  it  is  scheme  (7.8))  and  the  one  used  in  (7.7).  The  resulting  scheme  is  more 
accurate  than  scheme  (7.7)  but  less  robust,  we  have  then  (7.7)^  and 


(7.9)  i 


3  nn  +  1  0  r.n  _1_  1  „n-l 


2  nn  +  1  _  9  nn  -]_  i  nn_1  _n+I  _L.  i  +  l  i  -  n  +  1  yi_,n  +  l 

2  P  «•;  1  P<J+  2  pij  .  p.+lf  +  p.-li  +  Pi;  +  1  +  Pij— 1  -  4pi; 

f  -  -f- 


At 


!nn  +  1  4-lnn  +  1— nn+1  _i_lnn  +  1  I  n  n  +  1  1  ™n  +  1  r.  n  +  1  -L  1  r,n  + 1 

V+(M0.)-  2P ■'/ h  P-ii  +6P»-^  +  + 


lDn  +  1 +IDn  +  1_Dn  +  t  +lDn  +  l  ID"  +  1  +lDn  +  l_Dn  +  l  +1D"  +  1 

+  V+  (M,,) 


L=  (n+l)Al): 


we  shall  use  (7.8)  for  those  grid  points  at  distance  h  of  the  boundary. 
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8.  Further  Comments.  Conclusion 

Operator  splitting  methods  provide  a  systematic  way  to  obtain  efficient  solution  methods  for 
parabolic  equations  and  inequalities  in  high  dimension.  This  point  of  view  has  been  illustrated  by  the 
design  of  efficient  methods  for  fundamental  problems  such  as  the  Zakai  equation,  the  Navier-Stokes 
equations  and  time  dependent  obstacle  problems  (other  applications,  like  liquid  crystal  problems  are 
discussed  in  [8]).  Concerning  now  more  specifically  the  Zakai  equation  we  would  like  to  point  out  the 
following  facts: 

(a)  Due  to  the  very  high  space  dimension  associated  to  practical  problems,  there 
is  no  hope  that  fast  solution  can  be  achieved  for  several  years  to  come.  Indeed 
the  methods  described  here  are  quite  fast,  since  the  experiments  in  [8]  show 
that  the  solution  cost  per  time  step  is  of  the  order  of  the  number  of  grid  points; 
however  this  number  is  so  large  that  the  task  is  still  formidable,  requiring 
among  other  things  highly  massive  parallelism  and  huge  computer  memory. 

(b)  It  is  perplexing  to  realize  that  several  partial  differential  equation  specialists 
advocate  to  solve  deterministic  systems  of  advection-diffusion  type  (such 

as  the  Navier-Stokes  equations)  by  particle  methods  where  the  advection  is 
treated  via  the  integration  of  a  (large)  system  of  ordinary  differential 
equations,  and  the  diffusion  by  a  random  walk  method.  Indeed  the  stochastic 
dynamical  systems  leading  to  the  Zakai  equation  are  already  in  the  above  form 
since  they  couple  the  system  of  ordinary  differential  equations  modelling  the 
dynamic  of  the  problem  to  a  noise  ” generator”,  and  this  suggests  therefore 
to  deal  more  directly  with  the  original  form  of  the  problem. 

(c)  Since  the  Zakai  equation  is  linear  even  if  the  associated  dynamical  system  is 
not,  we  can  see  it  as  a  simplication.  This  is  a  quite  superficial  conclusion 
since  in  the  case  of  a  treatment  by  particle  or  characteristic  methods  we 
shall  have  nevertheless  to  solve  a  very  large  number  of  nonlinear  ordinary 
differential  systems  very  close  to  the  original  one  in  order  to  accurately 
construct  the  characteristic  curves  associated  to  the  advection  term. 


Our  very  conclusion  concerning  the  Zakai  equation  will  be  therefore  that  due  to  lack  of  computing 
power  it  is  hopeless  that  genuine  real  life  systems  can  be  treated  by  this  approach;  (at  least  for  several 
years  to  come)  in  fact  the  situation  is  very  close  to  the  one  encountered  in  the  solution  of  the  time 
dependent  Schrodinger  equation  where  only  the  simplest  molecules  can  be  investigated,  precisely  due  to 
the  dimension  difficulty. 
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